/*************************************************************************
 *                                                                       *
 * This file is part of the FreeMAPs distribution.  For full copyright   *
 * information, see COPYRIGHT and COPYING.LESSER.                        *
 *                                                                       *
 * Author:           sqbang                                              *
 * Email:            sqbang@ustc.edu                                     *
 * Date:             2013.12.10                                          *
 *                                                                       *
 *************************************************************************/

#ifndef FreeMAPs_include_COMMON_tbox_MathUtilities_h
#define FreeMAPs_include_COMMON_tbox_MathUtilities_h

#include "FreeMAPsConfig.h"

//C++ headers
#include <vector>
#include <limits>

//COMMON headers
#include "COMMON/tbox/Complex.h"

namespace FreeMAPs
{

/**
 * Class MathUtilities is a utility that provides some basic math-related
 * functions and routines for initializing data to signaling NaNs and
 * POSIX constants like INT_MAX, FLT_MAX, DBL_MAX, etc.  Signaling
 * NaNs force a trap if they are used in a numerical operation, so
 * they are a useful way to track uninitialized floating point data.
 * For example, setting integer values to INT_MAX is a useful way to track
 * uninitialized integer values.
 *
 * IMPORTANT: To properly trap operations based on signaling NaN values,
 *            the routine IEEE::setup_floating_point_exception_handlers()
 *            must be called.  This is normally done in the
 *            SAMRAIManager::startup() routine.
 *
 * The implementation of this class depends heavily on the particular
 * computer architecture and how it implements floating point arithmetic
 * and hardware traps.
 *
 * Note that the class @see tbox::IEEE also provides operations for
 * dealing with signaling NaNs.  This class provides the actual
 * implementation for such operations with functions in tbox::IEEE
 * calling the approriate operations provided here.  The class
 * tbox::IEEE is not templated on the data type and so calling the
 * operations provided there may be easier in some cases, such as in
 * codes built based on earlier versions of SAMRAI.
 */

template<class TYPE>
class MathUtilities
{
public:

    /**
     * @brief Return the value 0 for the template type.
     */
    static TYPE
    get_zero();

    /**
     * @brief Return the value 1 for the template type.
     */
    static TYPE
    get_one();

    /**
     * @brief Return the IEEE signaling NaN for the template type on
     * architectures that support it.
     *
     * Using this value in a numerical expression will
     * cause a program abort.
     *
     * For float, double, and Complex types, the usual signaling
     * Nan value will be returned if defined for the architecture
     * and compiler.  For char and int types, the POSIX max value for
     * the type is returned.  For other template types, zero is returned.
     */
    static TYPE
    get_signaling_NaN();

    /**
     * @brief Return true if the supplied value is NaN; else, false.
     *
     * For float and double, will check value against signaling NaN.
     * For Complex will return true if either real and imaginary part
     * is signaling NaN and false otherwise.
     * For other types will return false always.
     *
     * @param value Value to test
     */
    static bool
    is_NaN(
        const TYPE& value);

    /**
     * @brief Set vector entries to value given by get_signaling_NaN().
     *
     * @param vector SAMRAI array to set
     */
    static void
    set_vector_to_signaling_NaN(
        std::vector<TYPE>& vector);

    /**
     * @brief Set array entries to value given by get_signaling_NaN().
     *
     * @param array Pointer to first array entry
     * @param n     Integer to number of array entries.
     */
    static void
    set_array_to_signaling_NaN(
        TYPE* array,
        int n = 1);

    /**
     * @brief Return true if given values have a relative difference
     * smaller than sqrt(mach_eps) for the template type.
     *
     * Valid for float/double/Complex only.  For Complex will
     * return true if both real and imaginary part values
     * have a relative difference smaller than sqrt(mach_eps).
     *
     * For non float/double/Complex will just return (a == b).
     *
     * @param a
     * @param b
     */
    static bool
    equal_eps(
        const TYPE& a,
        const TYPE& b);

    /**
     * @brief Return max value for the template type.
     *
     * For boolean type, will return "true".  For Complex type,
     * will return a Complex value with both real and imaginary
     * parts set to the POSIX max value for type double.  For
     * other types, will return the POSIX max value for the type.
     */
    static TYPE
    get_max();

    /**
     * @brief Set vector entries to value given by get_max().
     *
     * @param vector SAMRAI array to set
     */
    static void
    set_vector_to_max(
        std::vector<TYPE>& vector);

    /**
     * @brief Set array entries to value given by get_max().
     *
     * @param array Pointer to first array entry
     * @param n     Integer to number of array entries.
     */
    static void
    set_array_to_max(
        TYPE* array,
        int n = 1);

    /**
     * @brief Return min value for the template type.
     *
     * For boolean type, will return "false".  For Complex type,
     * will return a Complex value with both real and imaginary
     * parts set to the POSIX min value for type double.  For
     * other types, will return the POSIX min value for the type.
     */
    static TYPE
    get_min();

    /**
     * @brief Set vector entries to value given by get_min().
     *
     * @param vector SAMRAI array to set
     */
    static void
    set_vector_to_min(
        std::vector<TYPE>& vector);

    /**
     * @brief Set array entries to value given by get_min().
     *
     * @param array Pointer to first array entry
     * @param n     Integer to number of array entries.
     */
    static void
    set_array_to_min(
        TYPE* array,
        int n = 1);

    /**
     * @brief Return epsilon value for the template type.
     *
     * For boolean type, will return "true".  For integer type, will
     * return "1".  For Complex type, will return a Complex value with
     * both real and imaginary parts set to the POSIX epsilon value for
     * type double.  For other types, will return the POSIX epsilon value
     * for the type.
     */
    static TYPE
    get_epsilon();

    /**
     * @brief Set vector entries to value given by get_epsilon().
     *
     * @param vector SAMRAI array to set
     */
    static void
    set_vector_to_epsilon(
        std::vector<TYPE>& vector);

    /**
     * @brief Set array entries to value given by get_epsilon().
     *
     * @param array Pointer to first array entry
     * @param n     Integer to number of array entries.
     */
    static void
    set_array_to_epsilon(
        TYPE* array,
        int n = 1);

    /**
     * @brief Return the minimum value of a and b.
     *
     * For Complex type will return the value with the minimum norm.
     *
     * @param a
     * @param b
     */
    static TYPE
    min(
        TYPE a,
        TYPE b);

    /**
     * @brief Return the maximum value of a and b.
     *
     * For Complex type will return the value with the maximum norm.
     *
     * @param a
     * @param b
     */
    static TYPE
    max(
        TYPE a,
        TYPE b);

    /**
     * @brief Return absolute value of a.
     *
     * For int, float, and double types will return the absolute
     * numerical value.  For other types, will return the input value.
     *
     * @param a
     */
    static TYPE
    abs(
        TYPE a);

    /**
     * @brief Return nearest integral value.
     *
     * If the argument is halfway between two integral
     * values then round away from zero for float and double types.
     * For other types, will return the input value.
     *
     * @param a
     */
    static TYPE
    round(
        TYPE a);

private:
    static TYPE s_zero;
    static TYPE s_one;
};

/*
 * Template specializations.
 */
template<>
bool
MathUtilities<float>::is_NaN(
    const float& value);

template<>
bool
MathUtilities<double>::is_NaN(
    const double& value);

template<>
bool
MathUtilities<Complex>::is_NaN(
    const Complex& value);

template<>
bool
MathUtilities<float>::equal_eps(
    const float& a,
    const float& b);

template<>
bool
MathUtilities<double>::equal_eps(
    const double& a,
    const double& b);

template<>
bool
MathUtilities<Complex>::equal_eps(
    const Complex& a,
    const Complex& b);

template<>
Complex
MathUtilities<Complex>::min(
    Complex a,
    Complex b);

template<>
Complex
MathUtilities<Complex>::max(
    Complex a,
    Complex b);

template<>
int
MathUtilities<int>::abs(
    int a);

template<>
float
MathUtilities<float>::abs(
    float a);

template<>
double
MathUtilities<double>::abs(
    double a);

template<>
float
MathUtilities<float>::round(
    float a);

template<>
double
MathUtilities<double>::round(
    double a);

} //namespace FreeMAPs

#include "COMMON/tbox/MathUtilities.inl"

#endif
